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Abstract 

Probabilistic cellular automata are prototypes of non equilibrium critical phe- 
nomena. This class of models includes among others the directed percolation 
problem (Domany Kinzel model) and the dynamical Ising model. The critical 
properties of these models are usually obtained by fine-tuning one or more 
control parameters, as for instance the temperature. We present a method 
for the parallel evolution of the model for all the values of the control param- 
eter, although its implementation is in general limited to a fixed number of 
values. This algorithm facilitates the sketching of phase diagrams and can 
be useful in deriving the critical properties of the model. Since the criticality 
here emerges from the asymptotic distribution of some quantities, without 
tuning of parameters, our method is a mapping from a probabilistic cellular 
automaton with critical behavior to a self organized critical model with the 
same critical properties. 

64.60.Ht, 64.60.Lx, 05.50.-hq 
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I. INTRODUCTION 



Recently, several papers have appeared discussing the relations between self or- 

ganized criticality (SOC) and usual critical phenomena. Some of them stress the 
fact that one can reformulate classic critical systems (namely, directed percolation) in a 
way indistinguishable from SOC, while others [||] focus on the role of the control and order 
parameters. 

We started our investigation from the observation |^ that one can express the problem 
of directed site and bond percolation in a form reminiscent of the invasion percolation 
process [^J or the Bak-Sneppen self organized model |TD|. The advantage of this formulation 
is that the critical value of the percolation probability does not need to be adjusted carefully, 
but instead emerges from the probability distribution of a set of continuous variables, while 
the original model is defined in terms of Boolean variables. 

The directed percolation problem can be formulated in terms of probabilistic cellular 
automata (PCA) [|ri|]. PCA are very general models that include for instance the kinetic 
Ising model fl^ . 

In this paper we show how any critical PCA may be mapped into a SOC model. The 
mapping is presented constructively in Sec. |I[ It can also be considered as a multi-site 
coding technique particularly adapted to probabilistic systems (where usual multi- 
site performs badly). From a computational point of view, this algorithm allows a quick 
determination of phase diagrams and computation of critical properties. In Sec. ^ we apply 
this method to the study of the Domany-Kinzel model of directed percolation and to the 
two dimensional Ising model. 

On the other hand, the mapping implies that to each PCA corresponds a SOC model 
defined in a high or infinite dimensional space. This correspondence can give some insight 
in the nature of the SOC phase, as addressed in Sec. ^V\. We end with some conclusions and 
perspectives. 



II. THE FRAGMENT METHOD 

We deal with probabilistic cellular automata, i.e., discrete models defined on a lattice. 
Let us consider explicitly the one dimensional Boolean case. A configuration at time t + 1 
is obtained from the configuration at time t by applying in parallel a probabilistic rule to 
each site. The rule is implemented on a computer by comparing (pseudo) random numbers 
with a certain number of fixed parameters (probabilities). One can think of PCA as the 
evolution of a deterministic discrete system on a random quenched field (the set of random 
numbers). 

For simplicity, we refer to the directed site percolation problem in 1-1-1 dimensions, where 
the higher p, the higher the probability of percolating. In this case one can visualize the 
random field as the height of a corrugated landscape, and p as the water level. There will 
be percolation if the water is able to percolate on the corrugated landscape, i.e., if the plane 
at height p is not completely blocked (in this directed model water is forbidden to back- 
percolate). For each value of p, we denote with a one the sites that are wet, and with a zero 
those that are dry. 



2 



One can stack a set of planes, and let them evolve in parallel. We can read the state of a 
certain site for all values of p as a vector of ones and zeroes, each component being labeled 
by p. 

In the initial configuration sites are either wet or dry, independently of p. Thus, all 
vectors are either filled with ones or with zeroes. Going on with the percolation process, a 
component p of the vector at a certain site and time t + 1 will be wet if there is at least a wet 
component at the same height among its neighbors at time t, and if the height of the random 
field at that position is less than p. One can easily express this in computer language. Each 
component corresponds to a bit in a computer word. With words of n bits, p can assume the 
values 0/n, 1/n, . . . , z/ra, . . . — V)/n. The (bitwise) OR of the words in the neighborhood 
gives a one for all values of p for which there is at least one wet neighbor. Given a random 
number r in that site, all planes with p > r have the possibility of percolating. This is 
expressed in computer terms by taking a word R{r) filled with zeroes up to a fraction r of 
bits and then with ones, and performing the AND of R{r) with the previous word. Iterating 
this procedure, we get in the last line of the lattice (say at time T) a set of partially filled 
words. If at time T a word has the bit number k equal to one, this means that for p = k/n, 
water would have percolated to that site (given the set of random numbers). 

The procedure can be generalized to words of arbitrary length. In the limit n —>■ oo, the 
Boolean vectors become the characteristic functions of subsets of the unit interval, that we 
call fragments. The manipulation of fragments is not limited to this bitwise implementation, 
as we shall see in the following. 

The fragment expressions do not depend explicitly on the control parameter p. The 
critical value of p and the critical scaling law of the order parameter are obtained a posteriori, 
from the distribution of fragments. 

Let us now formalize these concepts. For simplicity we refer to the Domany-Kinzel (DK) 



model |]TT|, which is a simple one dimensional PGA. We denote with x* = 0, 1 the state of a 
site i at time t, i = 0, . . . , L — 1, with L the size of the lattice. We shall simplify x' = x*"*"^, 
x± = xj^i- All space index operations are modulo L (periodic boundary conditions). The 
evolution rule may be written as 

x' = [r < p](x_ © x+) V [r < q]x^x+ (1) 

where © represents the eXclusive OR operation (sum modulus two), V the OR operation, and 
the multiplication (or A) stands for the AND operation, with the usual priority rules. The 
control parameters p and q are fixed, and r = r • is a random number uniformly distributed 
between zero and one, and where [logical expression] is one if logical expression is true and 



zero otherwise IIT^,p!5 



In the case of directed site percolation p = q, Eq. (|T]) can be rewritten as 

x'{p) = [r < p]{x^{p) V x+{p)). (2) 

where we emphasize the dependence of x on p. 

The fragment approach consists in reading xj{p) as the value of the characteristic function 
of the fragment X* at p. The expression [r < p] is the characteristic function of a fragment 
R{r) = [r,l). 

Equation (EI) in terms of fragments is 
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X' = i?(r)(X_ VX+) 



(3) 



that does not depend on p. The Boolean functions AND, OR, XOR, NOT correspond to 
the set operations intersection, union, symmetric difference and complement, respectively. 
We shall use the same symbol for Boolean and set operations. The initial configuration is 
independent of p, this means that are either the empty set or the unit interval, according 
with Applying the set operations we obtain the asymptotic fragments X°°. 

For a given value of p, is one if the point p belongs to the fragment X* and zero 

otherwise. Thus we can obtain the asymptotic value of xj{p) from the asymptotic fragments 
X*, that evolved without an explicit dependence on p. If some function of the x- exhibits 
a phase transition in correspondence of a critical value pc, this behavior can be extracted 
from the asymptotic fragments. For instance, the density p 

1 ^ 

p(p) = TE^r(p) (4) 
^ 1=1 

is proportional to the number of fragments the point p belongs to, 

p(p) = 7E[pe^n- (5) 

^ 1=1 

The above procedure can be applied to all probabilistic cellular automata. The practical 
recipe for the implementation is 

PCA: Express the model as a probabilistic cellular automaton whose evolution rule only 
uses Boolean expressions, and convert the control parameters {pi,p2, . . . ,Pm) to ex- 
pressions like [rk < Pk] where pk appear alone on the right side. 

Fragments: Replace the variables x- with fragments X* C [0, 1)™", and substitute [r^ < pk] 
with R{rk) {[rk > Pk] with its complement R{rk)). The initial configuration x° is 
replaced by X° = -R(x°). 

Implementation: Implement the fragments as arrays of bits (the simplest approach) or 
as sparse vectors (see later) and iterate the rule. 

Criticality: The asymptotic distribution of fragments gives the critical properties (control 
parameters and exponents) of the original model. 

Let us illustrate separately each of the previous points. 



A. PCA 



The evolution rule is generally expressed by means of transition probabilities. Note, 
however, that the transition probabilities do not completely characterize the problem for 
the damage spreading transitions [Q, since there are many ways of actually implementing 
the probabilistic choices in a computer code. The general approach for deriving a Boolean 
expression from transition probabilities is to write formally the future value of the dynamical 
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variable (the spin) as a function of the spins in the neighborhood and of the transition 



probabihties converted to random Boolean variables 



I.e., X 



fix- 



[r < j9ij, [r < 



.). Then there are several ways ||l8| , p!9| of expressing a Boolean function using a set of 
standard Boolean operations like AND, OR, XOR, NOT. Clearly, one should spend some 
effort in reducing the length of the resulting expression. Sometimes (see the Ising model 
in Sec. p|) one has to transform from [r < f{p)] to something like [/~^(r) < p\ (or more 
complex expressions). 



B. Fragments 

The method can be applied to any number of parameters. In the case of m parameters 
Pi, . . . the fragments are subsets of the m-dimensional unit hypercube. For instance, in 
the general DK model, Eq. (|l]), there are two control parameters (p and q) and the fragments 
are subset of the unit square. In this way it is possible to draw a sketch of a phase diagram, 
in just one simulation. 

However, if one is interested in crossing the critical surface along one line (see the compu- 



tation of critical exponents in Sec. [II D|) , one has to express the parameters as functions of a 
single variable, say s, and transform the expressions accordingly. For instance, the directed 
bond percolation problem corresponds to the curve q = p{2 —p), which can be expressed as 
p = s and q = s{2 — s) . The corresponding expression is 

x' =[r < s\{x- © x+) V [1 - a/1 - r < s\{x-X+) (6) 

which gives the following fragment expression 

X' = R{r){X^ © X+) V R{1 - VT^){X-X+). (7) 



C. Implementation 

The simplest way of implementing a fragment on a computer is by means of an array 
of n bits and using bitwise Boolean operations. Generally one uses computer words (32 
or 64 bits) for efficiency, but it is possible to use multiple words to increase the sampling 
frequency of probability. The numerical advantage over other multi-spin approaches is 
the use of just one random number for all the n simulations. Referring to a p-layer as a cut 
of the fragment space-time configuration with a given value of p, we see that the layers at 
different p are not independent, since they use the same random numbers. The influence of 
these correlations is discussed in Sec. [II D| . 

One can increase the sampling frequency around the region of interest (for instance the 
critical region) by appropriately deflning the correspondence of the bits with the values of 
p. This affects the way R{r) is implemented. With the fragment method it is still possible 
to perform simulations starting from a single site (Grassberger method) keeping track of 
nonzero fragments. The method is powerful if one uses a small interval around the critical 
point, so that all clusters for various p are similar. When using two parameters (say p and 
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q) one has to implement differently expressions like [r < p] (fill the unit square in the p 
direction from r to 1) from [r < g] (fill the unit square in the q direction from r to 1). 

The alternative approach in representing fragments consists in keeping track of the start- 
ing and ending points of all segments that form a one dimensional fragment. The approach 
is very similar to the treatment of sparse matrices, so we call it the sparse fragment method. 
The rules of combining sparse fragments are more complex than above. On the other hand 
in this way one has infinite precision, which comes useful in finding the critical behavior as 



explained in the Sec. IIP 



In general a fragment is formed by just one segment if the evolution rule can be expressed 
using only AND and OR, while for instance the XOR between two overlapping fragments 
causes holes. As an example, the site percolation rule Eq. @) can be implemented as sparse 
fragments by considering the evolution of the lower extremum a = a* of the segments [a, 1] 

as mm 

a' = max(r, min(a_, a+)). (8) 

Sometimes the problem can be reformulated without XOR operations. For instance, the 
bond percolation problem (^) can be rewritten as 

x' = ([r_ < V ([r_|_ < p]x+) (9) 

with two random numbers per site. The evolution of this rule can be easily implemented 
using sparse fragments. 



D. Criticality 

The critical properties of the original model are obtained from the asymptotic distribution 
of fragments. The fragment method introduces strong correlations among p-layers as also 
noted in Ref . P] . One can exploit these correlations considering differences in the p direction. 
If the patterns for different p-layers have similar sizes, the fluctuations cancel out. This 
happens in general if the rule does not contain XOR (see Fig. |l|). On the other hand, the 
XOR generally implies strong variations of clusters with p, so that the fluctuations can in 
principle be wider than uncorrelated simulations (see for instance Fig. |^). 

A powerful method for the computation of critical quantities exploits the scaling relation 

m(p, t) = a-^^^mia^/^ip - p^) + p^, at) (10) 

numerically solving it for the unknown (3, v and pc- This is an easy task for the sparse 
fragments approach, since one can obtain m(p, t) and m(p', od) with p' = a^^^ip — pc) + Pc 
for each value of v and Pc- For the bit approach one has to compute the value of the exponents 
and Pc so to make all data collapse on a single (smooth) curve. This can be performed near 
approximating the curves with polynomials (or any other fitting function), and minimizing 
the of the regression. 
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III. APPLICATIONS 



In this section we show some results of the fragment method apphed to classical problems: 
the determination of the phase diagram and critical properties of the Domany-Kinzel model, 
and the two-dimensional Ising model. 

The first example is the one- dimensional directed site percolation, i.e., the line q = p 
of the DK model, Eq. (^). A snapshot of part of the asymptotic fragment configuration is 
shown in Fig. |^, with the plot of the density p{p). As illustrated above, if one computes 
the XOR dilution (the line g = of the DK model, Eq. (JlD), the fragments decompose into 
several segments, as shown in the inset of Fig. 0. Correspondingly, the integrated density 
converges slowly to a smooth curve. 

The complete phase diagram of the DK model can be obtained in just one simulation 
by iterating two-dimensional fragments. The plot of the asymptotic density p{p, q) is shown 
in Fig. ^. It compares well with those obtained with other methods |TT] , p!B| , pT| . From the 
convergent behavior of the contour lines, the position [p = 1/2, q = 1) of the discontinuous 
transition for the density is clearly indicated. Near the corner (p = 1, g = 0) the surface 
becomes irregular: this is due to the prevalence of the XOR in eq. ([1|). One can also 
investigate the chaotic phase of the DK model by iterating two fragment configurations with 
the same random numbers. The Hamming distance between two replicas with a given value 
of p and q is the {p, g)-component of the density of the XOR between the two asymptotic 
fragment configuration. A plot of the resulting Hamming distance is shown in Fig. ^. Here 
one can notice a trace of the density phase boundary (near (p = 0.8, g = 0)), due to the 
critical slowing down. 

For the directed site percolation problem we found Pc = 0.7055(4), /5 = 0.21065(5), z/ = 
1.7195(5) for a system of size 10^, in the interval 0.7 < p < 0.71 for different values of 
a = 1024,2048,4096,8192. The agreement with previous measurements [§] is satisfactory. 
Moreover, we want to stress that these values were obtained with data coming from simula- 
tions of less than 20 minutes of CPU time on a 150 MHz PC running Linux |^ . 

As a second application, we consider now the kinetic version of an equilibrium system, 
the two dimensional Ising model with heat bath dynamics [^,^]. In Appendix A we show 
how to express the evolution equation of this model as a totalistic PCA, and how to translate 
its evolution in fragment language. In Fig. ^ we show the plot of the magnetization m(p) 
with respect to p = exp(— 2 J) for an Ising model with reduced interaction constant J. The 
transition is well characterized by plotting the second moment (standard deviation) of the 
magnetization as a function of p. We found Pc = .172 ± 0.002 and (3 = .11 ± .002, in good 
agreement with the exact values Pc = (-\/2 — 1)^ and /? = 1/8. 



IV. CRITICALITY AND SELF-ORGANIZED CRITICALITY 

We have shown how the DK model and the Ising model, which are PCA with critical 
behavior, may be mapped into fragment models with no control parameter, that is models 
that show self organized criticality (SOC). It is evident that the fragment method may be 
applied to any critical PCA. This result is summarized in the diagram of Fig. ^ The state 
xKp) may be obtained by evolving the PCA with a given p (labeled by fp in the diagram) 
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or by building the fragments Xf and evolving them with the fragment method {F in the 
diagram) that does not depend on p. Finally, by probing Xj with a p-layer, that is by 
checking if X* extends down to p, we recover xj. Although it is easier to think of one 
dimensional fragments, this result is valid for any number of control parameters. 

It is interesting to describe a "traditional" SOC model with the fragment language, 
trying to obtain the p-layer description that would make the SOC model correspond to a 
usual critical model. 



Let us discuss the one dimensional Bak-Snappen model |T0| with nearest neighbors in- 
teractions. In this model one starts from an array of real numbers Oj, i = 1, . . . , L uniformly 
distributed in the unit interval. One looks for the minimum of and replaces it and its 
nearest neighbors with newly generated random numbers, again uniformly distributed in 
the unit interval. The system auto-organizes so that the distribution of follows a power- 
law (with exponent 1), with a non-trivial avalanche distribution. We define now a parallel 
version of the previous model (which is not very efficient from a computational point of 
view). For the sake of simplicity, we divide the discussion in two parts: the research of the 
minimum and the actual evolution. 

Let us visualize the as the lower extremum of segments (fragments) Xj, and cut the 
configuration with a line at height p. The minimum is localized at site k for which there is 
only one intersection. It can be expressed using a Boolean variable 5*^, (a Kronecker 6) 



'5*. = V^*(p) A4(p) ' (11) 

assuming that the minimum is unique in the continuous p limit. 

The fragments Xi at and nearest to the minimum are replaced by segments of random 
length (-R(r)) 

= X* © (Al,^, V V Al,^,)iRirl) © X*) (12) 

where the fragment A*^^ is completely filled if 5*;, = 1 and completely empty if A* = 0. 
The evolution can be expressed on a j9-layer as 

xt^' =xl® (61,^, V 61, V 5*+i,,)(H <P]® a;*). (13) 

Similar but more complex expressions can be found for the invasion percolation process. 
One can see that in these "traditional" SOC models there are long range space interactions, 
and also interactions among p-layers (see Eq. (0)). We think that the second ingredient is 
the most important: if one knows how some quantity like the density varies with p it is not 
difficult to imagine a mechanism that automatically reaches the critical point. It is still to 
be proved that this is the actual mechanism of SOC. On the other hand, Eq. (^) shows that 
there exists space and p-local mechanisms that can be classified as SOC. 



V. CONCLUSIONS 

The fragment method can be considered both as a recipe for numerical studies of phase 
diagrams and as a mapping from criticality to self organized criticality. For what concerns the 
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first topic, the possibility of liaving a sketcli of the phase diagram without huge computation 
resources is useful in determining the position of the critical line. Numerical applications 
of the fragment method will be presented in future work. From the theoretical point of 
view, we think that the formalism presented in this work allows a clear characterization 
of the basic properties of self organized models, suggesting analogies between usual critical 
phenomena and self organized ones. 
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APPENDIX: THE ISING MODEL 

Let us start by considering the one dimensional Ising model. Its reduced Hamiltonian 
can be written as 

^W = -^li^^^m (Al) 



where cTj = 2xi — 1, and Xi = 0, 1. We choose the heat bath dynamics, P3| , ^ for which the 
probability r(x — y) of going from a configuration x to a configuration y that can differ 
from X in a certain number {ik} of sites is 

exp{-n{y)) 

^^"E'exp(-7^(y'))) ^^'^ 

where the sum in the denominator extends over all combinations of the differing sites 
This transition probability does not depend on x and satisfies the detailed balance principle. 

The configuration y is not limited to differ from x only at one site: the evolution can 
be applied in parallel changing all even (or odd) sites. Since the transition probabilities 
do not depend on the previous value of the cell, the space-time lattice decouples into two 
noninteracting sublattices: one with even sites at even times and odd site at odd times, and 
the complementary one. By considering only one sublattice, the neighborhood of the one 
dimensional Ising model is the same of the Domany-Kinzel model. The kinetic Ising model 
is just a totalistic cellular automaton (without adsorbing states). 

The local transition probabilities r(a;i_i,Xi+i Ui) can be computed from Eq. ( |A2| ) 
considering a difference in just one site. They are 

r(0,0^1)=p/(l+p) 
r(0,l^l) = l/2 
r(l, 0-^1) = 1/2 
r(l,l^l) = l/(l+p) 
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with p = exp(— 2J). 

It is convenient to introduce here the totahstic functions that take the value one if 
the sum of the variables in the neighborhood is k and zero otherwise. An efficent way of 
building these functions is described in [0]. For the Domany-Kinzel neighborhood they are 



Co = x_ V 
Ci = x_ © x+ 

The evolution equation for the Ising cellular automaton is 

x' =[r < p/{l+p)]co V [r < l/2]ci V [r < 1/(1 +p)]c2. (A3) 

Before using the fragment method one has to invert the [r < f{p)] {[r > f{p)]) expressions 
and substitute them with _R(/~^(r)) (i?(/~^(r))); p-independent expressions like [r < 1/2] 
trasform to R{0) or -R(O) according with r. We leave them in the equations with the 
assumption that a true value means -R(O) and a false value means -R(O). We finally find that 

X' = R (-^) Co+[r< l/2]Ci + R (^—^) C2. (A4) 



where now the Ck are fragments. 

For the 2D square Ising model with heat bath dynamics, the Hamiltonian is 

J 

^(x) = -- ^hj^i+hi + 
^ i,j=0 

with (Tjj = 2xij — 1. The lattice decouples again in two noninteracting sublattices. 

Repeating the procedure as above, one has again a totalistic cellular automaton. Using 
the totalistic functions Ck of four nearest neighbors, we find that 



X ' = i? U ^ Co V i? ( ; ) Ci V 



1' 

r < — 
2. 



Co 




yR ( C3 V i? I \r— I C4 (A6) 



where p = exp(— 2J). 



The Ck functions can be computed efficiently using the homogeneous polynomials Dj |19 



Co = Ci V C2 V C3 V C4 
C, = Di® 
C2 = D2® 
Cs = D, 
C4 = D4 



where 
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il)4 — X-^-^X^ X ^.X 
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FIGURES 



FIG. 1. The density p vs. the control parameter p for the directed site percolation problem, 
Eq. (|2|). The inset shows asnapshot of the first 40 segments X*, Eq. (^) after t time steps.. One 
simulation with L = 320 and t = 1000. The resolution is 480 bits. 

FIG. 2. The density p vs. the control parameter p for the XOR dilution, Eq. (|^) with g = 0. 
The inset shows a snapshot of the first 40 segments X* after t time steps.. One simulation with 
L = 320 and t = 1000. The resolution is 480 bits. Notice that the simulation reproduces well also 
the point p = 1, for which p = 0. 

FIG. 3. The contour plot of p{p,q) of the Domany-Kinzel Model, Eq. (|l]), for a lattice of 
L = 2000 sites and t = 4000. The resolution is 128 x 128 bits. White corresponds to p = and the 
contour lines are drawn at 0.1 intervals. 

FIG. 4. The contour plot of the asymptotic value of the distance between two replicas of the 
system evolving under the same realization of the noise. The parameters are those of the previous 
figure. 

FIG. 5. The magnetization m and the susceptibility Var(m) (fluctuation of the magnetization) 
for a 2-dimensional Ising model of size L = 100 x 100, averaging over 10 samples every 1000 time 
steps after a transient of 10000 time steps. The resolution is 32 bits. 

FIG. 6. The diagram showing the mapping from criticality to self organized criticality. 



13 



Bagnoli, Rechtman, Palmerini 



Algorithmic Mapping Critlcaiity Into Self Organized Critlcaiity 




figure 1 



Bagnoli, Rechtman, Palmerini 



Algorithmic N/lapping Criticality into Self Organized Criticality 



1.0 
0.9 
0.8 
0.7 
0.6 
p 0.5 
0.4 
0.3 
0.2 
0.1 
0.0 



^ ^ M ■■, I .1 "I ",■ I ■ ■■ ■■■I.,, : 

0.9 I'l .H ,'| ■I'lijjii.Jii":* |».!i-:t 

Q g ■ ■ ■ J ' !:M I'lV- . ■. I ■ ■ ■ 

0.7 

0.6 

X 0.5 

0.4 

0.3 

0.2 

- 0.1 

0.0 

10 20 30 40 



0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 

P 



figure 2 



Bagnoli, Palmerini, Rechtman Algorithmic Mapping Criticality 




Bagnoli, Palmerini, Rechtman 



Algorithmic Mapping Criticality 




figure 5 



Bagnoli, Palmerini, Rechtman 



Algorithmic Mapping Criticality into Self Organized Criticality 



X° = R{x°) 



X? 



fp 



F 



figure 6 



